function [HH,stdVec,flagNotPD]=...
estimaHessian(parMode,Y,funcModel,solveOpt,addSol,prpar,prss,stIn,gSteps)
% Compute the Hessian through Finite Point Difference
% Alejandro Justiniano Dec 28 2011
% Called from estimaOuter
% stIn must contain .LPosterior, Posterior at the value of ParMode 
if nargin < 9 || isempty(gSteps)
    gSteps=1e-5;
end
NP=size(parMode,1);
disp('Begin Hessian');
%% 1. Check Initial Density Matches Posterior at Mode 
fInitial=feval(@lmjPosterior,parMode,...
    stIn.parVec,stIn.parPosEst,funcModel,Y,stIn.trainVec,...
    prpar,prss,solveOpt,addSol,stIn.ssPosEst);
disp(sprintf('Starting Log Posterior for Hessian=%10.5f',fInitial)); 
if abs( fInitial - stIn.LPosterior ) > 1e-4
    warning('Starting Posterior for Hessian and Posterior at Mode do not coincide') 
end 
%% 2. Beging Hessian Routine 
HH=hessian(@lmjPosterior,parMode,gSteps,stIn.parVec,stIn.parPosEst,funcModel,Y,stIn.trainVec,...
    prpar,prss,solveOpt,addSol,stIn.ssPosEst);
disp('Done with Hessian')
try
    save Hessian
catch
    disp('Could not save Hessian')
end
%% 3. Attempt Inversion 
HH=reshape(HH,[NP NP]);
HH=generalized_cholesky(HH);
HH=HH\eye(NP);
[~,flagNotPD]=chol(HH);
if flagNotPD~=0
    warning('Hessian is not positive definite')
    flagNotPD=1;
end
stdVec=sqrt(diag(HH));